Adaptive Resolution Simulation of Liquid Water 
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We present a multiscale simulation of liquid water where a spatially adaptive molecular resolution 
procedure allows for changing on-the-fly from a coarse-grained to an all-atom representation. We 
show that this approach leads to the correct description of all essential thermodynamic and structural 
properties of liquid water. 
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Scientific problems in physics are usually classified by 
typical scales of measurable quantities such as energy, 
length or time. The appropriately determined scale to- 
gether with a tailored model and/or experiment is at the 
outset of a successful study. The specific properties under 
investigation then typically define the appropriate level 
of resolution. Using the optimal set of degrees of freedom 
(DOFs) guarantees efficiency, accuracy and avoids huge 
amounts of unnecessary detail, which might even obscure 
the underlying physics. This approach dates back to 
the very beginning of modern physics and finds its logi- 
cal continuation in systematic coarse-graining efforts for 
modern computational materials science and biophysics 
problems 0, 0, S 0]i where full blown all-atom simula- 
tions are often beyond the possibilities of current and 
near future computers. However, in many problems of 
materials science and biology different time- and length- 
scales are intrinsically interconnected, far beyond giving 
constant prefactors on the next coarser level of detail. 
Multiscale simulations are emerging as a promising tool 
for such problems SiBBi- Ideally, in addition to 
coarse-graining, one would like to make the step back to 
characteristic fine grained configurations as well [13, HH ■ 
Recently, an adaptive multiscale scheme (AdResS) has 
been proposed by some of us that also allows a free ex- 
change of particles between regions of different resolu- 
tions and to adjust locally the level of detail at will, 
while maintaining equilibrium with the fluctuating en- 
vironment [l^ . 

In this letter we address the above issue for liquid wa- 
ter, which continues to be a very active research field for 
the obvious reason, that all life exists in a water environ- 
ment. Although many classical models have been devel- 
oped for the simulation of liquid water 13, 14, 1^ 16| . 
there is still no unique model that can reproduce all its 
anomalous properties. Water is intrinsically multiscale, 
as it plays different roles at different scales. For instance, 
in the case of biomolecules in solution, at short length 
scales the hydrogen bonding with water governs the lo- 
cal shape and stability of folded biopolymers, whereas 
the "hydrophobic effect" drives the organization of the 



FIG. 1: (Color) On-the-fly interchange between the all-atom 
and coarse-grained water models, a) The explicit all-atom water 
molecule is represented at the left, and the coarse-grained molecule 
at the right. The middle hybrid molecule interpolates between 
the two, cf. text. A schematic representation of the full system is 
shown in b) , where a hybrid region connects the explicit and coarse- 
grained levels of description. All the results presented in the paper 
were obtained by performing NVT simulations using ESPResSo |19l | 
with a Langevin thermostat, with a friction constant F = 5ps~^ 
and a time step of 0.002ps at T^e/ = 300A' and p = 0.96(;/c?Tt^ (the 
density was obtained from an NPT simulation with Pre/ = 1 cutm). 
Periodic boundary conditions were applied in all directions. The 
box size is 94.5 A in the x and 22 A in the y and z direction. The 
width of the interface layer is 18.9 A in the x direction. 

system at longer time and length scales (as for instance 
the formation of membranes in aqueous solution). The 
approach presented in this letter represent a first step 
towards multiscale modeling of such complex scenarios. 
For fluctuating systems, such as adsorption of water on 
a hydrophilic surface, we envision a system in which 
an atomistic description is needed only near the sur- 
face, and a mesoscopic resolution can be used for wa- 
ter molecules further away from it, with molecules freely 
moving around without "feeling" the local level of reso- 
lution. 

To accomplish this goal, we first propose a new single- 
site water model that reproduces remarkably well the es- 
sential thermodynamic and structural features of water, 
as obtained by detailed all-atom simulations. In the sec- 
ond step we adapt the AdResS scheme to define a robust 
and physically accurate procedure to smoothly join the 
explicit atomistic and coarse-grained resolution regions. 
The resulting multiscale hybrid/mesoscopic model sys- 
tem 12|, [iTI, ll8| is composed of explicit and coarse-grained 
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molecules as presented in Figure [TJ 

Coarse-graining of water - Recent work has fo- 
cused on simplified coarse-grained models that can repro- 
duce qualitatively the all-atom center-of-mass (cm) radial 
distribution function (rdf) of water ^ EH, [ll, [21, [24 1 . 
To reproduce further structural properties, such as the 
orientational preferences for nearest neighbor configura- 
tions (i.e. for hydrogen bonding), previous one-site mod- 
els 2l|, |22| require a posteriori adjustment of the local 
environment |22l | . The one-site model of water presented 
in this letter can reproduce the structural and ther- 
modynamic properties of a widely used all-atom water 
model (namely, Rigid TIP3P [l3|) with remarkable accu- 
racy, without requiring any orientational corrections [25| . 
To construct the model we follow an iterative inverse 
statistical mechanics approach proposed by Lyubarstev 
et a l. liil (for alternative procedures see also Refs. [2| 
and [27| ) . The aim of this scheme is to numerically build 
a coarse-grained effective Hamiltonian that can mimic 
the behavior of a set of physical observables of the sys- 
tem under consideration (see i.e. Matysiak et al. 28|, [29| 
where this idea has been used to define a coarse-grained 
protein Hamiltonian "anchored" to experimental data). 
Additionally, in order to match the pressure of the coarse- 
grained to the all-atom model, after each iteration, a 
weak constant force is added to the effective force in such 
a way that the total effective force and potential energy 
are zero at the 7 A cutoff distance [13] ■ 

Figure [2] shows the final results obtained when this 
procedure is applied to a box of water molecules. A per- 
fect agreement between the all-atom and coarse-grained 
center-of-mass rdfs is reached with our optimized effec- 
tive potential (shown in the inset of the figure) after 8 
iterations. The effective potential has a first primary 
minimum at about 2.8 A corresponding to the first peak 
in the center-of-mass rdf. The slightly weaker and sig- 
nificantly broader minimum at 4.5 A corresponds to the 
second hydration shell. The combined effect of the two 
minima leads to a local packing close to that of the all- 
atom TIP3P water. Our effective coarse-grained poten- 
tial is quite different from the previously suggested po- 



tentials [20|, [2l|, [221: while in previous one-site models 
the deepest minimum corresponds to the second hydra- 
tion shell, the absolute minimum in our model is found 
in the first shell. In order to more thoroughly quan- 
tify the structural properties of our model (that are not 
completely defined by the rdf) we computed the angular 
distribution between the center-of-mass of three nearest 
neighbor molecules, and the distribution of the orien- 
tational order parameter q as defined by Errington et 
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where 



ipjk is the angle formed by the lines joining the oxy- 
gen atom of a given molecule and those of its nearest 
neighbors j and k. The parameter q measures the ex- 
tent to which a molecule and its four nearest neighbors 
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FIG. 2: (Color) a) The center-of- mass rdfs for explicit (ex) [cyan 
line] and coarse-grained (eg) [red line] regions of the hybrid system, are 
shown together with the rdfs corresponding to a bulk all-atom simula- 
tion (black line). The local O-H and H-H rdfs for the explicit molecules 
in the hybrid system compare equally well to the standard bulk simu- 
lations (not shown). The optimized effective potential for the coarse- 
grained model is shown in the inset as a function of inter-particle sep- 
aration, r. (b) The center-of-mass angular distribution between three 
nearest neighbors for the three cases studied as given in (a), (c) The 
analogous distributions of the orientational order parameter q. 

adopt a tetrahedral arrangement js^. The good agree- 
ment between the explicit and coarse-grained water mod- 
els shown in Figures[2](b) and (c) indicates that although 
our coarse-grained model is spherically symmetric and 
therefore does not have any explicit directionality, it ef- 
fectively captures the orientational preferences of hydro- 
gen bonding. 

As a consequence of the reduced number of DOFs, 
there is a time scale difference in the dynamics of the 
coarse-grained system, that is faster than what is pre- 
dicted from all-atom simulations. The diffusion coeffi- 
cient of the center-of-mass for the Rigid TIP3P model is 
D w (3.7 ± 0.3) X lO^^^i!. while for the coarse-grained 

model is £) « (8.1±0.4) x IQ-^^:^. While an accelerated 
time scale can be very advantageous in some cases, one 
can also adjust D by an increased background friction in 
the Langevin thermostat [3, [sij . 

Multiscale procedure - In order to smoothly join 
the atomistic and coarse-grained resolutions, we apply 
the AdResS scheme of Praprotnik et al. [l^, [l3| to cou- 
ple a model system composed of explicit and coarse- 
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grained molecules. Half of the simulation box is occupied 
by atomistic (Rigid TIP3P) water molecules while the 
other half is filled with the same number of correspond- 
ing coarse-grained molecules as schematically presented 
in Figure [1] The two regions freely exchange molecules 
through a transition regime, where the molecules change 
their resolution and their number of DOFs accordingly. 
The interface region contains hybrid molecules that are 
composed of an all-atom molecule with an additional 
massless center-of-mass particle serving as an interaction 
site. The transition, which needs to be smooth in order to 
be used in MD simulations, is governed by a weighting 
function w{x) G [0, 1] that interpolates the interaction 
forces between the two regimes, and assigns the identity 
of the particle. We used the weighting function defined 
in 12], in such a way that w = I corresponds to the 
atomistic region, and w = to the coarse-grained region, 
whereas the values < w < 1 correspond to the interface 
layer. The atomic and mesoscopic length scales are cou- 
pled via the intermolecular force acting between centers 
of mass of molecule a and P as: 



where FJ^*|™ 



T^atom 

.313 -icjfj 



9(7" 



is the 



J/3 dri^j^ 

sum of all pair intermolecular atom interactions between 
explicit atoms of the molecules a and /3 and F"'^ — 
— ^ is the corresponding effective intermolecular force 
between their centers of mass, r^^j^ ~ r^^ ^^j/i '^^ ^s*^" 
tor between atom i in molecule a and atom j in molecule 
/3 and Rq.^ = Rq — R/3 the vector between the centers 
of mass of molecules a and /3, with the corresponding X 
coordinates Xa and Xp. Each time a molecule crosses a 
boundary between the different regimes it gains or looses 
(depending on whether it leaves or enters the coarse- 
grained region) its equilibrated rotational DOFs while 
retaining its linear momentum. To supply or remove the 
latent heat caused by the switch of resolution this method 
is employed together with a Langevin thermostat . As 
discussed in an earlier publication, it is important to in- 
terpolate the forces and not the interaction potential if 
the Newton's Third Law is to be satisfied [18[. The key 
methodological issue in applying our approach to water 
is how to treat the long-range electrostatic interactions. 
We have chosen the reaction field (RF) method, in which 
all molecules with the charge center outside a spherical 
cavity of a molecular based cutoff radius Rc — 9A are 
treated as a dielectric continuum with a dielectric con- 
3, El. The Coulomb force acting on 



stant eRF [3J, |33|, 

a charge e^^ , at the center of the cutoff sphere, due to a 
charge e^^ within the cavity is: 



J/3 



.J0) 



1 



1 2{eRF-l) 



(2) 



There are three main reasons that make the RF method 
a natural choice to be applied with our scheme (Eq. [1]): 
the first reason is that it does not impose any artificial 
periodicity (that would interfere with the definition of 
the interface layer of the system) , the second one is the 
pairwise form of the reaction field term (Eq. [2]) , and the 
third one is that, as our scheme, the RF method works 
accurately only in a combination with a thermostat [s^l ■ 
The last point can be better understood if the switch- 
ing of the resolution is seen as a geometry induced phase 
transition Detailed comparisons between the bulk 

explicit simulations and the explicit regime in our hy- 
brid setup prove that this approach does not alter the 
structural properties of the water model studied. The 
width of the interface layer is determined by the maxi- 
mum range of electrostatic interactions, in order to allow 
the molecules to equilibrate on-the-fly while moving from 
one resolution to the other. In order to obtain a very 
smooth density distribution in the interface layer we use 
a conservative choice, that is, double the maximum range 
of interactions; however, a smaller interface layer can be 
used without affecting the fluxes between the layers, nor 
the results in the two resolution regimes. We observed 
that for a system containing only hybrid molecules, the 
pressure is increased in comparison to the reference all- 
atom system with the maximum deviation at w = 1/2 
(280% difference). This problem is solved by employing 
an interface pressure correction [l7l |. 

The structural (see Figure [2]) and thermodynamic 
properties of the bulk all-atom simulations are correctly 
reproduced in the hybrid system with the same mean 
temperature (0.1% difference) and pressure (0.5% differ- 
ence). The normalized density for the hybrid system is 
homogeneous in the coarse-grained and explicit regions 
with (very) small oscillations in the transition regime, cf. 
FigureO In order to prove the free exchange of molecules 
between the different regimes we have computed the 
time evolution of a diffusion profile for molecules that 
were initially localized at the interface layer. Figure [3] 
shows that these molecules spread out asymmetrically 
with time. This asymmetry arises from the aforemen- 
tioned difference in diffusion coefficient between the all- 
atom and coarse-grained regions. It is possible to obtain 
the same diffusional dynamics across the different reso- 
lutions by adding a local Langevin thermostat, as men- 
tioned above (36| . However, it is worth noticing that this 
difference in time scale can be advantageous for reaching 
longer simulation times in systems where multiple length 
and time scales are intrinsically present. 

In this letter a general and computationally efficient 
multiscale model for water has been presented. We 
have shown that a one-site coarse-grained model for wa- 
ter molecules can reproduce remarkably well thermody- 
namic and structural properties of the Rigid TIP3P wa- 
ter model, and that a smooth transition and free ex- 
change between coarse-grained and all-atom resolutions 
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FIG. 3: (Color) Top figure: normalized density profile in the x- 
direction of the hybrid system. Bottom figure: Time evolution of a 
diffusion profile for the molecules that are initially (at time t — ps) 
located in the interface region. The diffusion profile is averaged over 
fti 400 different time origins. Vertical lines denote the boundaries of 
the interface layer. 

is possible. We stress that the presented muhiscale ap- 
proach can be apphed to any other either flexible or 
rigid nonporalizable classical water model, e.g. SPC or 
SPC/E [35, 37, 38]. We envision that such a multiscale 
resolution of water will play an important role in the 
modeling of wet/dry interfaces and biomolecular simula- 
tions. 
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